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(54) Title: INDUSTRIAL PROCESS SURVEILLANCE SYSTEM 
(57) Abstract 

A system (10) and method for monitoring an industrial 
process and/or industrial data source (10). The system (10) includes 
a lime conrelation module (20), a training module (30), a system 
state estimation module (40) and a pattern recognition module (50). 
The system (10) generating time varying data sources, pixxessing 
the data to obtain time correlation of the data (20), determining 
the range of data, determining learned states of normal operation 
(30) and using these sutcs to generate expected values to identify 
a current state of the process closest to a learned, normal slate 
(40); generating a set of modeled data, and processing the modeled 
data to identify a date pattern and. generating an alarm (50) upon 
detecting a deviation from normalcy. 
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Industrial Process Surveiilajice System 



This invention was made with U.S. Government support under Contract No. 
W-3 1-109-ENG-38 awarded by the Depanment of Energy. The U.S. Government has 
certain rights in this invention. - ... 

The present invention is related generally to a method and system for carrying 
out surveillance of industrial processes using sensor or data source outputs. More 
particularly, the invention is concerned with a method and system for processing 
sensor data and using virtual data as an improved methodology over basic statistical 
approaches to industrial process surveillance. Further, the invention involves use of 
a plurality of techniques coupled for enhanced analysis of industrial process data. 

Conventional parameter-surveillance schemes are sensitive only to gross 
changes in the mean value of a process or to large steps or spikes that exceed some 
threshold limit check. These conventional methods suffer from either large numbers 
of false alarms (if thresholds arc set too close to normal operating levels) or a large 
number of missed (or delayed) alarms (if the thresholds are set too expansively). 
Moreover, most conventional methods cannot perceive the onset of a process 
disturbance, sensor deviation or data anomaly which gives rise to a signal below the 
threshold level for an alarm condition. Most methods also do not accoimt for the 
relationship between a measurement by one sensor relative to another sensor 
measurement. 

In another monitoring method, a conventional sequential probability ratio test 
("SPRT") technique has found wide application as a signal validation tool in the 
nuclear reactor industry. The SPRT method is a pattern recognition technique which 
processes the stochastic components associated with physical process variables and 
has high sensitivity for the onset of subtle disturbances in those variables. Two 
features of the conventional SPRT technique make it attractive for parameter 
surveillance and fault detection: (1) early annunciation of the onset of a disturbance 
in noisy process variables, and (2) the SPRT technique has user-specificablc false 
alarm and missed-alarm probabilities. SPRT techniques are primarily directed to the 
analysis of data from paired or multiple pairs of sensors in contrast to a large number 
of different process sensor data points. SPRT is also typically dependent on 
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assumptions of the data being independent of other dau sources and being Gaussian 
distributed data. The SPRT technique used alone therefore has certain shortcomings 
in identifying anomalies in processes. 

Other types of statistical techniques also have been developed for industrial 
process monitoring and analysis but have other insensitivities for certain classes of 
sensor data. 

It is, therefore, an object of the invention to provide an improved method and 
system for surveillance of industrial processes and apparati. 

It is another object of the invention to provide a novel method and system for 
on-line surveillance of industrial processes and apparati with multiple sensors. 

It is also an object of the invention to provide an improved method and system 
for evaluation of process data, on-line or off-line, from sensors or data accumulation 
sources. 

It is a funher object of the invention to provide a novel method and system 
for performing preliminary analysis of data for alarm conditions prior to data input 
to a SPRT system. 

It is an additional object of the invention to provide an improved method and 
system for masking selected sensor data and substituting virtual data to perform tests 
to determine whether abnormal process conditions or abnormal sensor conditions exist 
and whether or not to halt or modify the process under scrutiny. 

It is still another object of the invention to provide a novel method and system 
using training data characteristic of normal system and/or sensor and/or data source 
operation to compare with ongoing industrial processes and/or data accumulation. 

It is yet a further object of the invention to provide an improved method and 
system for processing data from a process to determine training data for normal 
operation, storing such training data on a computer storage media and analyzing real 
process data relative to the normal training data using a plurality of mathematical 
methodologies stored on a ROM or PROM storage medium. 

It is also an additional object of the invention to provide a novel method and 
system utilizing a virtual signal characteristic of normal state operation derived on the 
basis of correlation with a plurality of other process data values to compare with a 
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real process data signal set .for deriving the likelihood of an abnormal process or 
operation of data sources. 

It is yet another object of the invention to provide a novel method and 
apparatus lo accumulate training data to recognize any one of a plurality of specific 
states of operation and thereby identify a particular type of fault or condition present 
in a process or other system. 

It is also a further object of the invention to provide a novel method and 
apparatus for monitoring a process using training data to identify slowly changing 
operational sensor data characteristic of normal process changes. 

It is still an object of the invention to provide an improved method and system 
for determining whether a system or data source abnormality can be ignored without 
undesirable effects. 

Other advantages and features of the invention, together with the organization 
and manner of operation thereof, will become apparent from the following detailed 
description when taken in conjunction with the accompanying drawings described 
below. 

Brief Descri ption of the Drawings 
FIGURE 1 illustrates a schematic functional flow diagram of a preferred 

embodiment of the invention; 

FIGURE 2 illustrates a functional flow diagram of a time lead-lag conelation 

methodology; 

FIGURE 3 illustrates a functional flow diagram of a method of determining 
a full range of data by searching normal state training data; 

FIGURE 4 illustrates a functional flow diagram of a method for modeling 
behavior of commercial system operating slates; 

FIGURE 5 illustrates a functional flow diagram of a method for performing 
pattern recognition; 

FIGURE 6A illustrates sensor signals having a four second delay before 
applying a lead-lag method, and FIG. 6B illustrates the sensor signals after applying 
the lead-lag method; 

FIGURE 7A illustrates sensor signal data from pump 1 power with an SMSET 
estimate superimposed thereon; FIG. 7B illustrates the SMSET esrimation error 
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4 



between the SMSET estimate and the sensor signal data; and FIG. 7C illustrates a 
histogram of the error; 

FIGURE 8A illustrates sensor signal data from pump 2 power with an SMSET 
estimate superimposed thereon; FIG. 8B illustrates the SMSET estimation error 
between the SMSET estimate and the sensor signal data; and FIG. 8C illustrates a 
histogram of the error; 

FIGURE 9A illustrates sensor signal data from pump 1 speed with an SMSET 
estimate superimposed thereon; FIG. 9B illustrates the SMSET estimation error 
between the SMSET estimate and the sensor signal dau; and FIG. 9C illustrates a 
histogram of the error; 

FIGURE lOA illustrates sensor signal data from pump 2 speed with an 
SMSET estimate superimposed thereon; FIG. lOB illustrates the SMSET estimation 
error between the SMSET estimate and the sensor signal data; and FIG. IOC 
illustrates a histogram of the error; 

FIGURE llA illustrates sensor signal data for reactor outlet flow rate; FIG. 
IIB illustrates the SMSET estimation error between the SMSET estimate and the 
sensor signal data; and FIG. IIC illustrates a histogram of the error; 

FIGURE 12A illustrates sensor signal data for primary pump 2 flow rate; 
FIG. 12B illustrates the SMSET estimation error between the SMSET estimate and 
the sensor signal data; and FIG. 12C illustrates a histogram of the error; 

FIGURE 13 A illustrates sensor signal data for subassembly outlet temperature 
lAl; FIG. 13B illustrates the SMSET estimation error between the SMSET estimate 
and the sensor signal data; and FIG. 13C illustrates a histogram of the error; 

FIGURE 14A illustrates sensor signal data for subassembly outlet temperature 
2B1; FIG. 14B illustrates the SMSET estimation error between the SMSET estimate 
and the sensor signal data; and FIG. 14C illustrates a histogram of the error; 

FIGURE 15A illustrates sensor signal data for subassembly outlet temperature 
4E1; FIG, 15B illustrates the SMSET estimation error between the SMSET estimate 
and the sensor signal data; and FIG. 15C illustrates a histogram of the error; 

FIGURE 16A illustrates sensor signal data for subassembly outiet temperature 
4F1; FIG. 16B illustrates the SMSET estimation error between the SMSET estimate 
and the sensor signal data; and FIG. 16C illustrates a histogram of the error; 
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FIGURE 17A illustrates sensor signal data for reactor outlet temperature 
1534CF; FIG. I7B illustrates the SMSET estimation error between the SMSET 
estimate and the sensor signal data; and FIG. 17C illustrates a histogram of the error; 

FIGURE 18A illustrates sensor signal data for primary tank sodium level 530 
Float; FIG. 18B iilusiraies the SMSET estimation error between the SMSET estimate 
and the sensor signal data; and FIG. 18C illustrates a histogram of the error; 

FIGURE 19A illustrates sensor signal data for primary tank sodium level 531 
induction; FIG. 19B illustrates the SMSET estimation error between the SMSET 
estimate and the sensor signal data; and FIG. 19C illustrates a histogram of the error; 

FIGURE 20A illustrates standard deviation of SMSET errors for each of the 
data in FIG. 7-19; and FIG. 20B illustrates the mean value of SMSET errors for each 
of the data in FIG. 7-19; 

FIGURE 21 shows subassembly outlet temperature ("SOT") and SMSET 
estimates and in particular FIG. 21 A illustrates time dependent normal SOT for 3F1 
in the EBR-II nuclear reactor; HG. 21B illusuatcs normal SOT for 3C1; FIG. 21C 
illustrates normal SOT for 5C2 and FIG. 21D illustrates normal SOT for 7A3; 

FIGURE 22A-D illustrates SMSET estimation error for each of the data of 
FIGS. 21A-D, respectively; 

FIGURE 23A-D illustrates SPRT results for each of the data of FIGS, 21A-D, 
respectively; 

FIGURE 24 A corresponds exactly to FIG. 21 A; FIG. 24B includes a linear 
drift component compared to FIG. 21B; and FIGS. 24C and 24D correspond exactly 
to FIG. 21C and 21D, respectively; 

FIGURE 25A conesponds exactly to FIG. 22A; FIG. 25B includes the effect 
on SMET estimation error of the linear drift of FIG. 24B; and FIGS. 24C and 24D 
correspond exactly to FIGS. 22C and 22D, respectively; 

FIGURE 26A corresponds exactly to FIG. 23A; FIG. 26B illustrates the 
SPRT results for the linear drift error of HG. 24B; and HGS. 26C and D 
corresponds exactly to FIG. 23C and D, respectively; 

FIGURES 27A and 28B corresponds exacUy to FIGS. 21 A and 21B, 
respectively; FIG. 27C includes a temporary amplitude pulse of 0.25% of the signal 
magnitude; and FIG. 27D corresponds exactly to FIG. 21D; 
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FIGURES 28A and 28B corresponds exactly to FIGS. 22A and 22B; FIG. 
28C illustrates SMSET estimation error for the amplitude pulse effect of FIG. 27C 
and FIG. 27D corresponds exactly to FIG. 22D; 

FIGURES 29A and 29B corresponds exactly to FIGS. 23A and 23B; FIG. 
29C illustrates SPRT results of the amplitude pulse in FIG. 27C; and FIG. 29D 
corresponds exactly to FIG. 23D; 

FIGURE 30A illustrates EBRII subassembly temperature data 3F1 but includes 
a uniform gain change compared to FIG. 21 A and FIGS. 30B-D correspond exactly 
to HGS. 21B-D; 

FIGURE 31 A illustrates the SMSET estimation error for the gain change of 
FIG. 30A; and FIGS. 3IB-D correspond exactly to FIGS. 22B-D, respectively; and 

FIGURE 32 A illustrates the SPRT results for the gain change of FIG. 30A 
and SMSET analysis of FIG. 31A; and FIGS. 32B-D correspond exactly to FIGS. 
23B-D, respectively. 

Detailed Description of Prefe rred Embodiments 
The system 10 herein includes a methodology (see FIG. 1) and apparatus for 
surveillance of sensor or data accumulation configurations in industrial, utility, 
business, medical, investment and transportation applications. The system 10 is useful 
for sensitive identification of the onset of sensor or data source degradation, process 
or system anomalies, or the onset of change to a different operational state. The most 
preferred form of the system 10 comprises a synergistic integration of four techniques 
to provide enhanced surveillance capabilities compared to conventional approaches 
(including neural networks), and also provide enhanced reliability and improved 
computational efficiency. The four elements that make up the most preferred 
surveillance form of the system 10 are embodied in four different methodologies 
generally characterized as a time correlation module 20, a training module 30, a 
system state estimation module 40 and a pattern recognition module 50. 

Many attempts to apply conventional neural networks to nuclear, 
petrochemical or any industrial process surveillance applications have met with poor 
results in part because they fail to take into account lead-lag relationships (lack of 
proper time correlation of the data sets) between the various sensors or data sources. 
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In one example, a pipe is instnimented with a sequence of N thermocouples ("TCs") 
which measure the instantaneous temperature of the fluid passing through the pipe; 
and the signals from these TCs are displaced in lime along the fluid stream flow path. 
If the sampling inter\'al for the sensors is smaller than the transit time for fluid 
passing through the pipe, any attempt to monitor these signals with a neural net will 
produce very high uncertainties in the estimated signals for the TCs. For example, 
if it takes ten seconds for the fluid to pass through the pipe, and the sampling interval 
is one second, then at any given instant in time, TC(N), at the outlet end of the pipe, 
is seeing fluctuations that passed TC(1) ten seconds ago. These TCs may still contain 
a small degree of correlation due to gross changes in fluid temperature from a heat 
source or sink that is upstream of the pipe; however, the more valuable intersensor 
correlation that arises from local temperature penurbations carried along the pipe will 
be lost. This same phenomenon degrades the performance of neural networks and 
other pattern-recognition paradigms applied to any processes wherein the physical 
sensors or data sources are displaced in time across the process they are monitoring. 
Other examples of time delays in correlated systems include: systems with slow data 
flow rates and/or large, physical distances (oil refineries, power plants, HVAC 
systems, and financial systems), delays due to analog or digital electronics (noise 
fillers and large capacitors) or transmission delays (satellite communications, or 
transmitting data over different BUS systems. 

In a preferred form of the time correlation model 20, a Leadlag component 
of the invention (see FIG. 6) performs dynamic, real-time intersensor lead-lag 
adjustments. The Leadlag module 20 performs adjustments so that the output signals, 
which are then input subsequently into the SMSET routine (the system state 
estimation module 40), are optimally correlated and impart the maximum information 
content to the pattern recognition module 50. The Leadlag module 20 is attached 
hereto as a computer software Appendix A. The Leadlag module 20 accomplishes the 
adjustment function by performing, for each pair of signals, an iterative regression 
procedure that generates a vector of correlation coefficients with respect to lag time. 
This vector of correlation coefficients is a unimodal concave function of lag time. 
Thus, the optimal lag time between the pair of signals is identified simply by 
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searching for the zero-crossing of the first derivative of the vector with respect to the 
lag nmt. 

In other forms of the invention it may be unnecessary to utilize the Leadlag 
module 20, as noted in FIG. 1 wherein the option exists to skip use of the module 20. 
This could occur for systems in which there is adequate time correlation or if lime 
shifting is not needed to achieve correlation. If the Leadlag module 20 is not utilized 
or the data has already been processed by the Leadlag module 20, the data is 
preferably input to a training module 30. 

In a preferred embodiment this training module is a MiniMax module 30 
which searches through all the observations for all signals or data during a training 
lime period to construct training vectors that include the highest point and lowest 
point for each signal or data space under surveillance. A computer software Appendix 
B sets forth the MiniMax module 30. The MiniMax module 30 produces an 
"optimal" training set. It is optimal in the sense that it contains only, at most, 2N 
vectors, where N is the number of signals or data points in the system; and these 
vectors span the full range that all sensors or data sources have noted during the 
available training period. Wherever two or more sensors or data sources 
simultaneously attain maxima or minima, the resulting number of training vectors will 
be less than .2N. 

In another form of the invention both the Leadlag module 20 and the MiniMax 
module 30 can be skipped, and the data can be input directly to the system state 
module 40. 

Once the MiniMax module 30 has consu^cted a system model (or been 
skipped as noted above), the system state estimation module 40 (such as the preferred 
Subsystem Multivariate State Estimation Technique CSMSET") module) models the 
behavior of a system through examples of the operating states of the commercial 
system being modeled. A computer sofflvare Appendix C sets forth the SMSET 
module 40. In general, the system state estimation module 40 can be any one of a 
variety of modeling methodologies, such as auto regressive moving average, a neural 
network, or a Kalman filtering technique or an empirical methodology. 

The SMSET module 40 utilizes its memory of the learned states of the 
commercial system in conjunction with a single new observation to provide an 
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estimate of the current •true" system st2ie. States of the system are represented by 
vectors whose elements are comprised of direct values of system parameters 
(measured signals) as well as any transformation of these system parameters that 
produce scalar values, e.g. calculated parameters based upon measured data. The 
SMSET module 40 does not require the slate vector elements to be linearly 
independent as do most other types of estimation techniques. The learning process, 
which results in a "learned-state" matrix, is performed according to the MiniMax 
module 30 and the Leadlag module 20 described hereinbefore. 

The basic methodology of the SMSET module 40 involves the input of a new 
observation of the behavior of a system that is compared with the "memory" of 
previous system behavior embodied in the leamed-state matrix. A series of 
mathematical operations are performed that generates an estimate of the states in the 
system^s memory that is "closest" to the new observation. The definition of "closest" 
that is used by the SMSET module 40 is the slate that is lying closest to the new 
observation from the point of view of a set of rules that determine the association of 
two vectors. From this closest state, an estimate of the "true" state of the system is 
performed for each and every element of the state vector. Thus, given a set of 
current observed parameters of a system, the SMSET module 40 provides an estimate 
of the current true state of the system. The value of this method is that an estimate 
of all of the values of the system parameters in the state vector can be provided even 
if the current observation vector is incomplete (e. g. some sensors or data sources may 
have failed or are no longer available), contains erroneous or faulty elements (some 
sensors may have drifted, become uncalibrated, become contaminated with high noise 
levels, etc.), or even if the new system state does not coincide with previous 
operating states. However, the new system state must, in a general sense, be bounded 
by the domain of the states used to develop the system memory (leamed-state 
matrix). 

This estimation of the u^e current state of the commercial system, including 
estimated values of all system parameters, is used in conjunction with the actual 
measured system parameters to ascertain the opcrability of sensors (or other data 
sources) and disturbances in the system state. This state estimation process can 
further be described as* an inference engine that accepts as input a set of learned states 
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and a new observation of the commercial system. After a series of operations are 
performed by the inference engine on this input, the result is an estimate of the 
learned state "closest" to the new observation. The definition of "closest" used here 
is the state lying closest to the new observation from the point of view of a set of 
rules that determine the association (overlap) of any two vectors. Another result is 
the estimation of the "true" value of each and every element in the new observation 
vector in the form of an estimated state vector. The series of operations performed 
in the inference engine consist of various matrix operations. First, all pairs of learned 
states are preferably associated two at a time using a rule set to create the elements 
of a recognition matrix. Next, the new observation is associated with each learned 
state using the rule set to produce a vector that has the same number of elements as 
the number of learned states. The largest element value in this vector identifies the 
"closest" learned state to the new observation. Finally, the normal matrix product of 
this vector with the recognition matrix produces a set of linear combination 
coefficients for combining the learned states into the estimated state vector. This 
methodology, when applied to any true state of a commercial system that is a 
combination of the learned states, yields a very close approximation to the true state. 
The actual closeness achieved depends most upon nonlinearities arising from the rule 
set and physical and/or random fluctuations in the variables and is demonstrated by 
direct testing. General experience with use of this method for real operating 
commercial systems has indicated predictive capabilities typically no worse than 4-A 
0.5% ad normally +A0.1%. 

Once the SMSET module 40 has modeled the data as described hereinbefore, 
the data is input to a pattern recognition module 50, such as the Sequential 
Probability Ratio Test {"SPRT") module. The computer software for the SPRT 
module 50 is in Appendix D. This SPRT module 50 is a sensitive pattern recognition 
method that can detect the onset of subtle degradation in noisy signals with high 
reliability, and with quantitative false-alarm and missed-alarm probabilities. Output 
from the SMSET module 40 is provided as a set of estimated signals (also called 
"virtual signals") for each sensor under surveillance. These virtual signals arc fed into 
a network of interacting SPRT modules 50 together with the actual sensor readings. 
Each of the SPRT modules 50 receives one sensor-signal, vinual-signal pair. If any 
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sensor degradation or process disturbance starts to affect the output of one or more 
signals under surveillance, the SPRT module(s) 50 provide an annunciation to the 
operator and an actuator signal to the control system, which can selectively as needed 
automatically swap in the virtual signal to replace the degrading sensor signal, or data 
source. Further details of the SPRT module 50 are described in USPN 5.459.675, 
which is incorporated by reference herein. 

The above-described combination of methodologies enables identification of 
a faulted process, a particular type of fault, a faulted sensor or data source or faulty 
data itself and enables actions to be taken to correct or modify the process being 
monitored. 

In some cases when a failed sensor, or improper data stream source, is not 
important to the continued operation of a commercial system, the user can continue 
operating the commercial system or process if the sensor or data source were 
operating normally. For example, the system 10 can operate to substitute in a 
modeled estimate into a actual commercial system or process as input to replace a 
failed sensor or failed data source. This allows the commercial system or process to 
keep operating. 

Since the system 10 does not rely on analytical modeling by itself, it is 
applicable to a wide variety of processes and systems, such as petro-chemical, power 
generation, automotive, manufacturing, medical, aeronautical, financial and any 
system in which signals are available for processing that are related to the 
commercial system process operation or performance. The only requirement of the 
system 10 is that there is some type of cross-correlation, be it linear or nonlinear, 
between the signals used as input to the system 10. The signals can be linear, 
nonlinear, stationary, nonstationary, clean or noisy (with an arbitrary distribution). 
The system 10 uses a database of historical operation data to model the commercial 
system or process. The database is assumed to contain data from all relevant 
operating modes of the system; however, if a new mode of operation is encountered 
and is determined not to be a result of commercial system or sensor failures, a new 
vector can be added to the existing training matrix to incorporate the unanticipated 
operating mode in the system model. 
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The following nonlimiiing examples illustrate various aspects of the invention 
described herein. The data used is all taken from the EBR-II reactor at Argonne 
National Laboratory (West). 



The effect of time delay correlation is described in this example, and FIGS. 
6A and 6B show two voltage signals with a four second delay between them. Before 
the signals are processed by the Leadlag module 20 (see FIG. 6A), the correlation 
coefficient is 0.0182 which implies no correlation versus processing through the 
Leadlag module 20 to obtain a conelation of 0.9209 (sec FIG. 63). When the set of 
signals, or data, being used is more than two, all the possible pair combinations are 
used to calculate maximum. possible correlation coefficients so all signals can be 
properly correlated. 



An experiment to determine the accuracy of the invention (the "SMSET" 
methodology generally) was carried out using sensor data from the Experimental 
Breeder Reactor II (EBR-II) at Argonne National L^oratory (U558). The sensor data 
set contained 13 signals from sensors monitoring EBR-II. Table I shows the SMSET 
Estimation accuracy for EBR-II Data. Table I includes the channel numbers and 
descriptions for each of the sensor signals used in the analysis. The experiment was 
conducted in three steps; first the SMSET module was trained using two days worth 
of EBR-II data, next the trained SMSET module was used to estimate the state of 
approximately 110 hours worth of EBR-II data, and then the accuracy of the 
estimates was analyzed. For each of the sensor signals listed in Table I, FIGS 7-19, 
respectively, show the sensor signal (top plot) and SMSET estimate superimposed, 
the middle plot shows error between the SMSET and the sensor signal (in percent of 
the signal magnitude), ad a histogram (bottom plot) of the error. The histogram plots 
are compared to a Gaussian distribution with the same mean ad variance to give an 



Example I 



Examplg n 
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idea of how Gaussian the error signals. FIG. 20 provide a summary of the data of 
HGS. 7-19. 

A methodology entitled MiniMax (Appendix B) was used to train the system 
using the two days of training data cited above. After the MiniMax method was 
applied, a training matrix was constructed consisting of twenty-five unique vectors 
constituting an empirical model of the overall system. After creating the model, the 
methodology was then applied to the signals listed in the accuracy table. Each signal 
in the system has its own estimation error signal that is a measure of how close the 
pattern recognition model is representing the system 

relative to the sensor readings. The second column of Table I lists the standard 
deviation of the estimate error for all of the signals in the experiment in terms of 
each of the signals' magnitude. The magnitude of the signal is defined by its mean 
during normal operation. The third column in Table I lists the mean of the estimate 
error for all of the signals also in terms of the signal magninide. In general the 
estimate error standard deviations are in the range of 0.01 % to 0. 1 % and the esumate 
error means are centered around 0. Bar graphs of the tabular information are shown 
in FIGS. 20A and 20B as graphic representation of the accuracy information. 

Tablg I 



SMSET Estimation Accuracy Table for EBRII Data: 



Channel Number and 
Description 


SMET Estimate 
Error Standard 
Diviation (% of 
Sensor Magnitude) 


SMSET Estimate 
Error Mean Value (% 
of Sensor Magnitude) 


1) Primary Pump (f\ Power 


0.05245 


0.01241 


2) Primary Pump #2 Power 


0.14352 


0.06595. 


3) Primary Pump #1 Speed 


0.01078 


- 0.00001 


4) Primary Pump ffl Speed 


0.01272 


-0.00278 


5) Reactor Outlet Flowrate 


0.09585 


0.04452 


6) Primary Pump #2 
Flowrate 


0.06034 


-0.02495 


7) Subassembly Outlet 
Temperature 2B1 


0.04635 


0.00339 
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8) Subassembly Outlet 
Temperature 2B1 


0.04904 


-0.00118 ' 


9) Subassembly Outlet 
Temperature 4E1 


0.05664 


-O.UUJOo 


in\ Ciihacc^mhlv Outlet 

Temperature 4F1 


0.04926 


-0.00413 


11) Reactor Outlet 
Temperature 1534CF 


0.04727 


0.00513 


12) Primary Tank Sodium 


0.02440 


-0.00280 


12) Primary Tank Sodium 
Level 531 Induction 


0.00615 


0.00316 



Example III 

In FIGS. 21-32 examples of different sensor failure modes are shown along 
with how the system reacts to the failures. The preferred method of FIG. 1 is applied 
to the data. The sensor signals used in these examples are from a subset of 22 sensor 
signals used in the system. The 22 sensors monitored the EBR-II subassembly system 
at Argonne NaUonal Laboratory (West). Each of FIGS. 21-32 contains four subplots 
in which the upper most plot is related to Subassembly Outlet temperature ("SOT") 
3F1, the upper middle plot is related to SOT 3CU the lower middle plot is related 
to SOT 5C2, and the bottom plot is related to SOT 7A3. The system applied in each 
of the examples uses the same training matrix, which consists of 83 vectors selected 
from a training data base containing almost a weeks worth of data taken once every 
minute. 

In FIGS. 21-23 are shown the results of using the system 10 during 
approximately 5.5 days of normal operation of EBR-II. FIG. 21 shows the SOT 
signals with their corresponding SMSET estimates (signal being the circles and the 
lines being the estimate). FIG. 22 shows the respective raw estimate errors (not in 
terms of the signal magnitude) derived by taking the difference between the SOR 
signals and corresponding SMSET estimates. Finally in RG. 23 the results are shown 
from applying the decision making module of the system 10 (the SPRT module 50-- 
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see Appendix D) lo ihe SMSET estimation errors of FIG. 22. The SPRT plots show 
a total of only three false alarms which is a false alann rate of 9.4 x la*. and this 
is weU within the specified false alarm rate of 1.0 x 10*\ 

One type of failure mode that is common among sensors is a slow drift. This 
type of failure can be difficult to detect early on especially when the sensor signals 
contain noise. FIGS. 24-26 illustrate a comparative example of processing data from 
this type of failure and failure identification. Signal #2 (FIG. 24B) has a 0.2% linear 
drift in it's mean over the 2.75 day period starting at 4000 minutes into the signal. 
The other sensors are operating normally. FIG. 25 shows the resulting SMSET 
estimation errors for each sensor signal. The error plot for signal W (FIG. 25B) 
shows evidence of drifting after the sensor signal has drifted approximately 0.05%. 
In FIG. 26 the SPRT method has determined that #2 (FIG. 26B) is drifting after 
approximately 0.05% of drift and that all other sensors are operating normally. 

Another type of failure that can occur is a step change in the sensor signal. 
This can be the result of a short in the sensor or DAS, a calibration error or for a 
variety of other reasons. FIGS. 27-29, show a example of this type of failure for the 
SOT measurements. In this example sensor signal #3 (FIG. 27C) contains a pulse 
with an amplitude of 0.25% of the signal magnitude. The pulse starts at 4000 minutes 
and lasts for 2000 minutes. FIG. 27 shows the sensor signals and the SMSET 
estimates for the four SOT signals. FIG. 28 shows the resulting SMSET estimation 
errors. The error signal for ff3 (FIG. 28C) shows that there is a problem starting ai 
4000 minutes and ending at 6000 minutes. The error signals are fed through the 
SPRT module 50, and the results are plotted in FIG. 29. Clearly, there has been a 
disturbance in sensor #3 (FIG. 29C) beginning at time 4000 minutes and ending at 
6000 minutes. 

In FIGS. 30-32 an example of a failure mode related to the sensor gain is 
shown. In this example the gain of the sensor signal changes over tune, i.e., the 
amplitude is increasing over lime. The gain begins changing linearly over time from 
a beginning value of 1 to a final value of 1+0.075% of the sensor magnitude. The 
system 10 for the estimation error is applied to the signals, and the results are shown 
in FIG. 31. A human operator would most likely not be able to tell that there is a 
problem even after 8000 minutes by looking at the sensor signal. In FIG. 31 A, it is 
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apparent that signal #1 is operating abnormally. This is confirmed in FIG. 32A by the 
SPRT results, showing a steadily increasing number of SPRT alarms over the 8000 
minute period. 

While preferred embodiments of the invention have been shown and 
described, it will be apparent to those skilled in the art that various changes and 
modifications can be made without departing from the invention in its broader aspects 
as set forth in the claims provided hereinafter. 
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Appendix A 

Computer software for Leadlag module which performs dynamic, real-time 
iniersensor lead-lag lime correlation adjustments. 
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/♦ IXADLAG OPTIMIZATION MODULE 



Irinciude <stdio.h> 
ItxncXudtt <mach.h> 
«include <stdlib.h> 
trincXude cstring.h> 

«de£iae pi 3.1415926 
#de£ine N 1024 
#de£ine xcsize 2047 
#de£ine fil^xcsize 2079 
»de£xne filsize 32 

void phAse_Bhift opcimizaciondnt argc.char argrvCl 120], inc method); 

void conv(float £(] , float g[]. float c[]. int size) ; 

void der(float yC), float diffll. int size) ; 

void roots(float di£fxc(], int size, float Toot) ; 

void fliplr (float a[], int size) ; 

void pm_ir*f o (char 



void mamdnt argc, char ••argv) 

{ 

int i, methods- 
char tempargvdO] 120) ; 

if(argc ;-5) { 

pm_inf o ( argv 1 0] ) ; 
exittl) ; 



£orli«0 ; i<argc;i++) 
atrcpy Itempargvli] .argv(il ) ; 

printf ("Enter l or 2 below. (\*1\" means the employment cf the") ; 
printfC derivative technique to find the shift, while \-2\''") ; 
printf (" means the application of direct max. correlation-); 
printf (" technique. ): Nn") ; 
scanf ("Vd" . dmethod) ; 

phaie^shif t^optimization (argc. tempargv, method) ; 



void phase shift ootimization(int argc, char argv[] [20). int method 
{ 

int i.j,ph; 
ifloat root; 

float filxcaizel .f:(xcsize) .blfil^xcsize) ; 
float xc If il^xcsize) . 5 il_xc (f il_xcsize) ; 
float difffil xc(fil_xcsize-l) ; 

FILE •ir.filelT*in£ile2.*infile3.*outfile. •outfilel. •outfile2: 
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AMENDED CLAIMS 

[received by the International Bureau 24 October 1997 (24.10,97); 
original claim 1 amended; new claims 17-30 added; 
remaining claims unchanged (5 pages)] 

1. A meihod for monitoring at least one of an industrial process and 
industrial sensors, comprising the steps of: 

generating time varying data from a plurality of industrial sensors; 

processing the time varying data to effectuate optimum lime correlation 
of the data accumulated from the plurality of industrial sensors; 

searching the lime correlated data to identify maximum and minimum 
values for the data, thereby determining a full range of values for the data from the 

industrial process; 

determining learned states of a normal operational condition of the 
industrial process and using the learned states to generate expected values of the 
operating industrial process; 

comparing the expected values to current actual values of the industrial 
process to identify a current state of the industrial process closest to one of the 
learned states and generating a set of modeled data; 

processing the modeled data to identify a pattern for the data and upon 
detecting a deviation from a pattern characteristics of normal operation, an alarm is 
generated. 

2. The method as defined in Claim 1 wherein the industrial process 
comprises a physical process. 

3. The method as defined in Claim I wherein the industrial process 
comprises a financial process. 

4. The method as defined in Qaim 1 wherein the step of dete rmining 
optimum time correlation comprises comparing pairs of sensor signals, each 
characteristic of a separate sensor and calculating a cross correlation vector over 
lime, applying a low pass filter to remove noise from the cross correlation vector and 

. determining phase shift between the sensor signals. 
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17. A method for monitoring at least one of an industrial process and an 
individual date source, comprising the steps of: 

sensing time varying data from at least one industrial data source of 

an industrial process; 

determining learned stales of a desired operational condition of the 
industrial process to use the learned states to generate expected values of the 

industrial process; 

comparing the expected values to current sensed values of the industrial 
process to identify a current slate of the industrial process closest to one of the 
learned states and generating data characteristic of the current state; and 

processing the data that is characteristic of the current state to identify 
a pattern for the data and upon detecting a deviation from a pattern characteristic of 
the desired operational condition, a signal is generated indicating at least one of the 
industrial process and the industrial data source is not of the desired operational 
condition. 

18. The method as defined in Claim 17 further including the step of 
searching the time varying data, before comparing the current actual values to the 
expected values, to identify minimum and maximum values for the data, thereby 
establishing a full range of values for the data. 
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19. The method as defined in Claim 18 wherein each said industrial data 
source is characterized by two data values associated with the minimum and 
maximum values. 

20. The method as defined in Claim 17 wherein said step of processing the 
data characteristic of the current state to identify a pattern comprises applying a 
sequential probability ratio test. 

21. The method as defined in claim 20 wherein data that is characteristic 
of the current state is processed to generate a set of modeled data which is further 
processed to identify the pattern for the data. 

22. The method as defined in Claim 17 wherein the industrial process is 
selected from the group consisting of a manufacturing process, a physical process, 
a chemical process, a biological process, an electronic process and a financial 
process. 

23. The method as defined in Claim 17 further including the step of 
substimting an estimated signal for said industrial data source upon deteaing the 
deviation from a pattern characteristic of the desired operational condition, thereby 
replacing a faulted data source enabling continued operation and monitoring of the 
industrial process. 
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27. The method as defined in Claim 26 further including the step of 
substituting an estimated signal for said industrial data source upon detecting the 
deviation from a pattern and characteristic of the desired operational condition, 
thereby replacing a faulted data source enabling continued operation and monitoring 
of the industrial process. 

28. The method as defined in Claim 26 further including the step of 
processing the time varying data to effectuate optimum time correlation. 

29. The method as defined in Claim 26 wherein the step of identifying a 
pattern for the time varying data comprises applying a sequential probability ratio 
test. 

30. The method as defined in Claim 26 wherein the industrial process is 
selected from the group consisting of a manufacturing process a physical process, a 
chemical process, a financial process, an electronic process and a biological process. 
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11. The method as defined in Claim 10 wherein the step of identifying a 
current sute of the industrial process closest to the learned state includes forming a 
combination of the learned states to identify a true state of the industrial process. 

12. The method as defined in Claim 10 further including the step of 
substituting an expected value for incomplete observations. 

13. The method as defined in Claim 10 wherein the industrial data source 
is selected from the group consisting of an industrial manufacturing process, a utility 
operation, a business operation, an investment process, weather forecasting and a 
transportation system. 

14. The method as defmcd in Claim 10 wherein the plurality of industrial 
data sources comprises a plurality of sensor pairs. 

15. The method as defined in Claim 10 wherein the step of processing the 
modeled data comprises applying a SPRT process. 

16. The method as. defined in Claim 10 further including a step of 
determining time phase shift between the plurality of the time varying data being 
output. 
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5. The method as defined in Claim 4 wherein the step of determining 
phase shift comprises differentiating the cross correlation vector with respect to lag 
time between each pair of signals and performing an interpolation to compute the root 
of the differential of the cross correlation vector. 

6. The method as defined in Claim 1 wherein the step of identifying a 
current sute of the industrial process closest to the learned state includes forming a 
combination of the learned states to identity a true state of the industrial process. 

7. The method as defined in Claim 6 further including the step of 
substituting an expected value for incomplete observations of the industrial process. 

8. The method as defined in Claim 1 wherein the step of detecting a 
deviation from normal operation comprises applying a sequential probability ratio test 
to the modeled dau. 

9. The method as defined in Claim 1 wherein the step of detecting a 
deviation from normal operation comprises performing a pattern recognition analysis 
using computer means. 

10. A method for monitoring at least one of an industrial process and an 
indusuial data source, comprising the steps of: 

generating time varying data from a plurality of industrial data sources; 

determining learned states of a normal operational condition of the 
industrial process to use the learned states to generate expected values of the 
operating industrial process; 

comparing the expected values to current values of the industrial 
process to identifying a current state of the industrial process closest to one of the 
learned states and generating a set of modeled data; 

processing the modeled data to identify a pattern for the data and upon 
detecting a deviation from a pattern characteristic of normal operation, an alarm is 
generated. 
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^JJY^^] Claimed is: 

1 . A method for monitoring an industrial process and industrial sensor, 

comprising the steps of: 

generating time varying data from a plurality of industrial sensors; 

processing the time varying data to effectuate optimum time correlation 
of the data accumulated from the plurality of industrial sensors; 

searching the time correlated data to identify maximum and minimum 
values for the data, thereby determining a full range of values for the data from the 

industrial process; 

determining learned states of a normal operational condition of the 
industrial process and using the learned slates to generate expected values of the 
operating industrial process; 

comparing the expected values to cunent actual values of the industrial 
process to identify a current state of the industrial process closest to one of the 
learned states ad generating a set of modeled data; 

processing the modeled data to identify a pattern for the data ad upon 
detecting a deviation from a pattern characteristic of normal operation, a alarm is 
generated. 

2. - The method as defined in Claim 1 wherein the industrial process 
comprises a physical process. 

3. The method as defined in Claim 1 wherein the industrial process 
comprises a financial process. 

4. The method as defined in Claim 1 wherein the step of determining 
optimum lime correlation comprises comparing pairs of sensor signals, each 
characteristic of a separate sensor and calculating a cross correlation vector over 
time, applying a low pass filter to remove noise from the cross correlation vector and 
determining phase shift between the sensor signals. 
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char cii. bu£tX61; 

>0; £lag-Oi i*0) jaOy 



vhiXft ( (c&i-'\a' ) b&(chi«roy) ) ( 

aaa-(t (iMp»c«(ch)});/-tKchi-- .MtAichl-'-O) 
/• t£ Cch— (k**; aoviD-j; j»Oi> •/ 
if (BOB— 1){ 

bof [i] -ehi 
bttf li4ii-'\0'j 

A«pl« Ij] • (float) strtod(bu£«fiIZLL) ; 

l*Oi 

£Ug«0: 

> 

) 

/• £clo»«(lo); •/ 

rfttxm (0)i 

) 



Katio (float ami. float nurnl) { 

if Kmal — 0) II ( B»&3 — 0 1) 

itLuxu O.O; 
alaa if (sn&l >• ma2) 

racom snal/Btalr 

•lam 

ra Lmu wmL/aamlt 

\ AtfZaSO 

iat ma« flag, it j« 

char c^. buf[lCly 

mx8»0i filagvQt i*Oi j*0} k»Oi 

vMla ( (eh«fg«te (Tnia) ) I •BOT) { 

wmm ( 1 (itapaea Ce&) ) b& (ehU' . ' ) && (ehl-' *M } i 
if (oh»ii'\A') {k^^i &o«sO-j; j«Oi} 
if (ina»l){ 

bttf til meht 

flag«ls 

} 

if ((m»— flUt(nag—l)){ . 



wo 97/49011 



PCTAJS97/10430 



41 

it { (Spo»l>-M) I 1 (SMgl>«BB) ) tubl-l? 
if ( [Spoft2>«IB) I I (S&sg2>*BB) ) iub2>l: 



•printf tsprtposl. •H£\ii- . Spoil) i 
fprdatS (tprrpoaa^ •\t\ti* , Spota) s 
fpriatf (»prta«gl# •^t\a* * Simgl) ? 
fprlnrf (vprtii6g2. •%f\it"*Sna92) } 

fseanf (f Xovl. "%f \zi* . btnpl) f 
faciaf {£««tl, \&' . &t«ap2} i 

If (lubll 

Cprlatf (cambl, •^£\n» , t«»p2) ; 

fprintf (ccadbl, "^fXa" . t»apl) i . 

fiunf CfXov2.*%t\n".tiC«Bpl) ; 
fiacuf C£«at2,«%f\&*«fctaap2} ; 

12 («^) 

fprlntf lctsab2,-\f\a*.tBap2) ; 

fprlatf (coB±i2,'%f\n«»t«^l) ; 
l£tf3oi4((donbl») (cat+1-^1) r (doubl«)p»r) •-0) ( 
gat03cr(47»ll) ; 

cpriatf (-n-a-af V ((float J <cat*i*l) )/( (float) L«agcliJ •100/2) ; 

} 

} 

gatox7(47.11) i 

cpriatf(»S-3.3f •.100.0) j 
gocoacy (34.15) : 
t«xrX»ooJc(VSZrZ«UB) I 

cprijitf (• air MT xzr •) ; 
tfotoiy (45»13) ; 
gatehd 1 

fclcsAivprrpofl) ; f clot* (tprtnagl) i fdott (fasti) : f cloaa iarrarl) i 

f aloaa (enbl) i f oloaa (£lo»l> ; 
fcloaa(Bprtpoa2)) £cloaa<»prtMg2) , fdoaa (faat2) ; fclo»*(arror2) ; 

f oloaa (CGBib2) ; f cloaa (flov2) i 

f oloaa (la) i 



raadixputSaopla (Till 
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fg«tt (faaf ^13.«zTor2) ; 

} 

£i»«kl«rrerl,0LrO} i 

for (i.»0ii<c«lcli«ngt^ii»»){ 

Cgats (btt£»13.«rroxl) ; 
tap»4toi (bttf ) I 

^vlucAl^* (tBsp-Manl) *ftvqi*Maal) / (olcLangtli'Il : 
f g«ta (taf , 13 • •zxor2) ; 

) 

/• Actul $nr c&lcul&tion ftad signal svsppiag logic */ 
«iihl»0» ■ub2»0i 

gl «X1 / vmxl&ncft 1 1 
}O-FXC70S*flqrt(vttrdazicaa) i 

xjk-log (irrx/ (1 -u^psx) > i 

BB-log( (1-Br=l) /JLIiPSU ) 

If (1—0) fprl&rf (pKrmms,*«if\a-,A^} 

1X(1*-13 £prdjitl(pmxws*«^f\&«,BB} t 

If (1—3) fprlatf (p«r»ma,-%f\n«.JU) I 

1£{1— 31 fprlac£(paramj,«\£\B«,BB)i 

) 

fcloa« <p«ms) ; 

Sposl*0; Spot3»0r Sn«gl*Oi Snag2«0i 

for (l«OpicL«ngxb:i*«) { 

/* fg«ta(tuf .13»«rror3) ) 
t«ap*Aco£ (bu£) : 

•/ 

fsevL£(ftrxorl.«%£\n«.ftt«Bpl) ; 
f tcuLf (arxar2. "%£\n- « fctci^)2) i 

KSpotl — JU)||(Spotl SB)) Spotl-Oi 
if ((Spos3 iJL)||(Spoa3 B8)) Spo«2-0i 
If {(BMgl — AA) 1 1 (Ssagl SB) ) Siugl»Or 
if ( (SnAg3 JU) 1 1 (fiugl •» BS) ) 5a*g3-0 1 

•tapposl- •gl*(Xl/2 • tupl) : ttopnmgl* *gl*(Kl/3 «- s«Bpl) ; 
i%«ppot3« -g2*(XB/2 - t8B^) ; fr«paog3« •g2*aC2/2 ♦ tWBp2) ; 

SpoBl *« atttpposl; Sa^gl *-« tttpMgli 
Spoa3 fCcppofl2i S&ag3 trspnagai* 

i£(9posl c. AX) Spoil«AA: if (Sa«gl <- XA) Sasgl-AAi 
if (Spes3 <- AA) Spo«3-AAr if(Saftg3 <m lA) Sbo93«AA: 

if (Sped >- BS) Sposl»BB: iKSiugi >. bb) Snft9l«aai 
if (Spoa2 >- BB) Spo»3«B&i i£iB&A92 >• &B1 Sa«92«BBi 
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Appendix D 

Computer software for performing pattern recognition by delecting onset 
degradation in noisy signals. 
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H function [s, nortt_di»t) - vprprod (x, y. •Ipha, bet*); 
V 

\ Thl» function performs the nonlinear operetion on aatricoo x and y at tha 

% heart of the K5CT. It ueet the following nonlinear 

S operator to coopsre vectors in the satricea; 
% £(di8t) - 1 / (1 4. alpha dlst^ata) . 

% where dlit ij tho noraalixed diatance between a vector in y and each 

% obsanration vector in x. 

\ 

\ If the conatant alpha and beta are not epecified, default values of 1 ara 

\ used for both of tham. 

\ 

\ Thia version of vprprod ratuma the noro of tho diatance between the 

% each vector in y and the ajtaioplar vetora In x. The vector distancee are 

% norsalised with raapact to the aua of the distance vectora. 

%. 

S 

s 

function tz. aorsa^diat) - vprprod (x. y. alpha, beta); 

if (nargin 2) 
alpha • Ij 
beta - 1: 

and 

[Ux.Xx] - aiae (x) ; 
tMy*MyJ - aiaa(y); 



if (MX — My) 

fprintf (' \n\n in vprprod (x.y) , the nuaber of coluana in x Buat\n' ) ; 

fpriatf (' equal the number of rows in yVn\n*); 

atr* ( ' coluane in x: msD2atrOCc), ... 

' , and the nuaber of rowa in y: nuDL2str (Hyl ] ; 

•rxor(atr); 

end 

\ Initialise arrays, 
z - caroa (Ux.My) t 
DOTS m xeroa (Ux.My) ; 
yd o xeroa (Hx*Mx) ; 
diat m xeroB(l7x,l)} 
tsp > zeros (Ux,Nx) ; 

for i • ltl*y 



% Rxpand each column of y into a Kx by Ky {- Hx by Mx) matrix. 
- I../. j\ * — 
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\ Resc^ V matrix to input matrix, 
y » Input ( : , 2 :size (Input, 2) ) ' ; 

% Calculate error. 

Err o Y (prt_colB, : ) - X_hat; 

% Calculate rms error for variables specified by prt^cols. 

Erms • diag(Err ♦ Err* ) ; 

Erms « sgrt(ErtnB' ./ (N - M - i) ) ; 

% Add timestamp to X^hat and Err arrays and transpose X_hat and Err arrays. 
X_hat «= [Time' ; X_hAt) ' ; 
Err « [Time' ; ErrT' ; 

if (print_flag >■ i) 
btm-fix (clock) ; 

f printf ( * \n\n Total time to complete MSET using %d training vectorsNn' , ... 
size(D«2)); 

f printf (' to analyze \d observations with %d variables each: N, M) ; 
fprintf ( ' %6 . 2f mins\n\n' , toc/60) ; 
end 
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XJbAt - xeroKilength(prc_cols) .K) 

cop > x«roa (M.l) ; 
Ermi - teroi(leagrh(prt_col8) .1) ; 
Err . 2«rosl length tprr_coXs) .N) ; 
DtY - teroBiH-lengthtouc^cols) .1) ; 
W ■ zaroB(K.l) ; 

% Main loop, calculate estim.ce vectors (X.hat) tor each input state 

% vector Y. 
for i-l:N 

DtY - vprprod(Dj)roj ' . Y^roj ( : , i) . alpha. b«t*» ; 
W . DDI • DtY; 

m U J >um(W) ; 
ii-£lxxd(W<thrsh) ; 
if (leftgth(ii)>0) 

¥{ii)»2erofl<lttngth(ii) .1) ; 
\ M m W./sum(V) ; 

tmp ■ D • W; 

Xhat(:,i) - tmp(prt_coi») : 

% Print me»«*ge if print_flag is let. 
if (print_fl»g >- i) 

if t (remix. floor(K/100)l 0) | ii--!)) 

fprintf ('\b\b\b') ; 

fprintf ('%3.0f ' , i/H»100» ; 

end 

end 

% scale the X hat matrix from the range 0:1 back to the range specified by 
% the rangeL and rangeU vectors, if neccaaary. 
if llangthlrangaL) — M) t (langth(rangeU) M) 

% scale each obaorvation vector in the training matrix. 

%RHN-lcngth(rangeUlprt_col«) ) ; 

%rU.rangeU(prt_cola) ; rL-rangeLlprt_cola) ; 

%X hat- (rut :) •ones 11, K) - . . . 

% " rL(:)*oneB(l.N)l ... 

% • X hat + rL( :)*anaa {1»N) J 

for i . l:lengtbTprt_coli) . i•^^^ ^ 

Xhatli.:) . (rangeU(prt_coiati)) - rangeLlprt.col. (i) ) ) ... 
X_hat(i.:) ♦ r»ngeLtprt_cola(i) ) ; 

end 

% Shift the X_h*t matrix b.cJc co th. ruige ot the data by addlag raageL to 
% the data, if necessary. 

if llengthlrangeL) — M) k (rangeU 0) 

% Shift each observation vector in the training matruc. 
tQIH-length(prr_cola) ; 

X^hat-X^hat ♦ rangeLlprt^cols) •ones (l^WW) ; 
%ior 1 ■ l:length(prt_coXs) 
. \ X^hatli.:) -X_hAt(i.:) * rangcL(prt_cols(i)); 
%end 

end 
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if nargxn alpha- l; beta-i; p_flag-0; end 

if nargin 8. beta.l; p_flag-0; end 

if nargin B, p.flag«0; end 

if prc^cois 0, prc_colfi-l :M; end 

if alpha -» 0, alpha«l: end 

if beta «• 0, beca-1; end 

\ Set print flag for debug printout. 
princ_flag « p_flag; 

% Print meflBage if print^flag is set. 
if (print_flag >» 1) 

fprintt ( ' \n Starting MSET \n\n' ) ; 

end 

% Scale the input matrix to Che range 0:1. if the rangeL and rangeU 
V vectors are specified. 

if (length (rangeL) M) k ( length trangeD) M) 

I Scale each observation vector in the training matrix, 
y. (Y-rangcL* ♦ones 11, K) )./... 

(rangcU' *onea tl.N) -rangeL' ♦ones (l.N) ) ; 

%for 1 - l:M 

% Yd,:) - (Y(i,:) - rangeL(il) / (rangeU(i) - rangeL(i)); 
lend 

V Print message if print_flag is set., 
if (prLat_flag >• 1) 

fprintf { • Data is scaled to 0:1 range .... \n\n' ) ; 

end 

end 

% Shift the training matrix by the values specified in rangeL. if the rangeL 
% vector is specified and rangeU is 0. 
if (lengtht rangeL) M) L (rangeU 0) 

* Shift each observation vector in the training matrix. 
Y - Y - rangeL* •ones (l.N) ; 

%for 1 . 1:M 

\ Yd, :) - Y(i, :) - rangeL(i) ; 
lend 

\ Print message if print^flag is set. 
if (print^flag >- 1) 

fprintf {' Data is shifted by the values in rangeL. ... \n\n' ) ; 

end 

end 

\ It the are output columns specified, then remove these variables from the 

% D and Y matrices. 

Djro3 • D; 

Y_pr03 • Y; 

if (out^cols — 0) 

Ojroj (out_cols. :) • [] ; 

YjJroj (out^cols, :) ■ (); 

end 

\ Print message if print^flag is set. 
if tprint_flag >• 1) 

fprintf:' percent complete: 0' ) ; 

end 
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\ tunccion lErma, X^hat, ErT»WW2.WWl - m»ct (Input. rangeL. rangeU, D. DOi. 



\ ouc cols. prt_cols. alpha, bata. p_f lag. thrah) ; 
% 

\ This function estimaccs che acate of a syacem using the (MSET) . It 

I allowB scaling o£ the input data to the range 

\ 0 - 1 and the uae of the algorithm m a predictive mode. 

% 

^ Function outputs: 

k Exms • Root mean squared errors for variables specified by 

I prt_col8 . 

* X_hat » Estimated states for variables specified by prt_col8. 

% Err - Estimation error for variables specified by prt^cols. 
% 

% Function inputs: 

% Input - input state matrix: an N by M4l array with the first colxian 

I containing a timestamp. The array contains n observation 

% vectors, with m dapandent variables in each vector. 

% rangeL « Length M vector, where each element sets the lower limit of 

% data for a corresponding variable. 

\ rangeU = Length M vector, where each element sets the upper limit of 

S data for a corresponding variable. 

\ D ■ Training set. an M by N array. 

\ DDI - Inverse of the similarity matrix (from auto_trainnsa) . 

% out cols - optional vector specifying which columns represent output 

% variables predicted by the algorithm. These variables are 

\ not included in the evaluation of W. 

\ prt_cols - Optional vector specifying which columns to keep in the 

I output . 

% alpha » Multiplicative factor in the vprprod nonlinear operator. 

% beta - Power factor in the vprprod nonlinear operator. 

% P_fl*g ■ Print/don't print run progress data for 0/1. 

\ thrsh - Cutoff threshold for weighting vector w (default 0) 



% The rangeL anri rangeU vectors are used to specify the anticipated range of 
% the variables. Data is linearly scaled from the ranges specified by these 
\ vectors to the range 0:1. If the raxigeL vector specifies the lo%fcr limit of 
\ the data while rangeU is 0. then the code will shift the data by subtracting 
\ the rangeL values from each element in a corresponding column of the data. 
\ It both o£ the rangeL and rangeU vectors are set to a scalar value of 0, 
\ the data will not be scaled. 

function (Erias, X_hat. Err.WW2,wwl » maet (Input, rangeL, rangetJ. D, DDi. ... 

out_cols, prt_colB. alpha, beta, p_flag, thrsh) ; 
if (nargin c ll) , thrah-O.O; end 
tic 

\ Case in which not enough inputs are specified, 
if nargin < 5 

error {' Required arguments are: Input. rangeL, rangeU, D and DDi\n\nM ; 

and 

% Remove time stamp from input matrix and transpose input matrix. 

Time - Input ( : . 1) ; 

Y - Input (: .3;size(Input.2) )' ; 

(M. H] - aize(y); 

% If out^cols. prt_cols. alpha, beta, and p_flag are not specified, pick 
\ default values for them. 

if nargin 5. out_cols»0; prt_cols»'l :K: alpha-1: beta-i; p_flag-0; end 
if nargin o. prt_coic«l :M; alpha* l; beta»l; p_flag«0; end 
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Appendix C 

Computer software for modeling behavior for examples of the operating states 
of the commercial system. 
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nozmrideo { ) : 
cextLooktYTXLOW.BLUE) ; 



1 



void textLoo)c(int txt. int bg) ( 

x«3tcbac)cgrouxxd<bg) ; 
c«xccolor (txt) ; 

1 

void NQwKame(int DcpoB. int ypos, int whcl 

{ . - - - 



int 1. c 
char tcsqp [ 4 ] ; 

ttXtLook ( LI CH T GR AY . BIACK) ; 
90Coxy(xpo8-12,ypos) ; 
printff •); 
gotoxy(xpos-12,ypOB) ; 
i«Of 

maxi*-!; 

wtdle ((c-g«tche())!-13){ 
if (c«»8) (i-i-1; 

if (i<01 i«0; 

) 

else( 

maxisiiiax( i.maxi) ; 
iwitch (wfac) ( 

casa 1: training [ i ] cc ; braak; 
caaa 2: input [il"c; braak; 
caaa 3: tflS9[i]»c; braak; 
caae 4: tenpiij «c;braa)c; 
case 5: tesplij 'c.-braaJc- 
dafault: break; 

} 
) 



} 

switch (vt^c) C 

caae 
caae 
caae 



1: traiaing[maxi+ll«* \0* ; break; 
2: input fmaxitl] >' \0 ' ; break; 
3: templinaxi+l] = ' \0* ; 

FACTOR«( float) acof(texnp); break; 
4: teinp[xnaxitl1 s ' \0 ' ; 

Fl*atoi{tenip} ; break; 
5: 

temp(maxi+ll"'\0' ; 
F3 -atoi i temp ) ; break ; 
default: break; 



caae 



caae 



Fl-Fl-1; 
F2-F2-1; 
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) 

if (ch-=72)l 

yy«yy2; 
if tyy<ll) 

yy-21: 

1 



if ( (ch-«8»tt(yy!"2X) ) ( 

if (yy«.ll) 

NewNamelxx.yy, 1) : /* l«craining 

if (yy««il3) 

N«wName(xx,yy,2) ; /* 2=input 
if (yy»-15) 

N«wNairve(xx,yy, 3) ; /* 3=output 
if lyy"17) 

NawNajne<xx,yy, 4) ; 
if tyy««19) 

KtfwMame txx,yy« 5) ; 
yy«yy+2; 

) 

if (yy21) 
xx-42; 

else 

xx«53; 



if { (yyB=21)tt{ch=«131 ) stop»l; 
/• TESTING CONTENTS */ 

/* 

90C0xy (1,22) ; 

cprintf ( • ■ ) ; 

gotoxy (1.22) ; 

cprintf ( •training - %s" , training) ; 
gotoxy (1.23) ; 

cprintf (• •)? 
gotoxy (1.23) ; 

cprintf (■ input = %s\n" . input) ; 
gotoxy(1.2C) ; 

cprintf ( • • ) ; 

gotoxy (1.24) ; 

cprintf ( 'output = %8\n" . output) ; 
V 

gotoxy (XX, yy) ; 



void InitialScrettn{ ) 

{ 

c«xtbac}cground(CyAN) ; 
clrscr ( ) ; 

textLoo)c{RCD.BIACK) ; 
bighvideo ( } ; 
gotoxy (24 ,2) ; 

cprintf ( " • ) 

gotoxy (24. 3 ) ; 

cprintf (• KSET SIGNAL VALIDATION SYSTEM •) 
gotoxy (24 . 4 ) ; 

cprintf ( • • ) 
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void GatSettings t ) 



{ 



int ch. yy* xx. 
char 'bailf; 



stop; 



/• default Fl and T2 

Fl»l; 

r2-2; 



columns xor flow 1 and flow 2*/ 



t«xtLoo)c( YELLOW. BLUE) : 
c?r!ntf('*The default input and output file names are listed below. 

c?rintf (''uU %c and %c to move to selection. <bac)cspace> to change name. 

24.25); 

gotoxy(22,ll) ; 

textLook ( BLACK . LICHTGRAY ) ; 
cprintfC 1) training data :•) ; 

gotoxy(40.11) ; 
tcxtLoo)c(LIGHTCFAY. BLACK) ; 

cprintfC %a -.training); /♦ use gotoxy (51, 11) V 

gotoxy(22.13) ; 

textLoo)c (BLACK. LIGHTGRAY) ; 

cprantfC 21 input data :"); 

gotoxy(40,13) ; 

textLooJc(LIGHTGKAY. BLACK); 

cprintfC %■ '.input); /' use gotoxy(51. 131 . * / 

gotoxy(22.15) ; 

textLoolc (BLACK. LIGHTGRAY) ; 

cprintfC 3) STM :*); 
gotoxy(40.151 ; 
textLook:(LIGHTGRAY. BLACK) ; 

cprintfC 2 /* gotoxy (51. 151 V 

got03cy(22.17).; 

ttxtLook( BLACK, LIGHTGRAY) ; 

cprintf C 4 ) Column flow 1 : " ) ; 

gotoxy(40.17) ; 

textLook(LIGHTGRAY. BLACK) ; ^ 

cprintf C 1 ' * ^ 

textLooic (BLACK. LIGHTGRAY) ; 

gotoxy(22.19) ; 

cprintfC 5) Column flow 2C}: 

gotoxy(40.19) ; 

textLoo)c (LIGHTGRAY. BLACK) ; 

cprintfC 2 N; 
gotoxy(37,21) ; 

textLoolc (RED, LIGHTGRAY) ; , , 

cprintfC done '); /'use goto (42. 21) •/ 

/♦ make selection */ 

gotoxy (53 . 11) ; 

yyll; 

xx-53 ; 

stop"0; 

while (stopt*l){ 
ch«getch( ) ; 

if (ch»-80){ 

yy«yy*2; 
if (yy>21) 

yyll; 
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1-0: 

whiles (c=gecc(fl) ) !=*\n' ) { 



) 

•/ 

int ChooseTrainingt I 

{ 

int ch. yy, xx. stop; 
char 'buff; 

textLook ( YELLOW , BLUE ) ; 
cprintfr There are four different ways of extracting training data from 



gotoxy(8,B); . 
cprintfC the training file, 
gotoxy(28, B) ; 
textLook (13. BLUE) ; 
cprintf ( •%-12s • . training) ; 
9Otoxyi40. 8) ; 
textLook (YELLOW. BLUE) ; 
cprintf choose one below: 

gotoxy(30,ll) ; 

textLook (BLACK. LIGHTGRAY) ; 

cprintf (' 1) Use File as is '): 

gotoxy (30.13) ; 

cprintf!" 2] MinMax method "!: 
gotoxy(30, 15) ; 
cprintf {• 3) Equally Spaced 
gotoxy(3 0,17) ; 

cprintf!" 4) Specify states '); 

/♦ make selection */ 

gotoxy(48, 11) ; 

yy-ll; 

xx«48; 

StOp'O ; 

while (stop! >1) ( 
ch«getch( ) ; 

if (ch«-80){ 

yyyy+2; 
if (yy>17) 

yy»ll; 



gotoxy(l,22) ; 
cprintf (■ 
gotoxydf 22) ; 

cprintf ("Choice = %d* , ( (yy-9) /2) ) ; 
•/ 

gotoxy(xx,yy) ; 



•); 



(ch"72) { 
yy«yy-2 ; 
if (yy<ll) 

yy»17; 



if 



(ch«el3) stop»l: 

TESTING CONTENTS 



return l (yy9) /2) ; 
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cprintfC Extending training matrix 

extendD ( ex } : 

) 

void .UserSpeci f ied ( ) ( 
InitialScroant I ; 

?^r!^fr'ThU feature is not available m the DEMO version 

?^rin^f r'Jsing MinMax (hit a key to continue) 

getch( ) ; 
HisMAXt ) ; 

) 

void BxtendDlfloat 6x» { 

Soat'm(MAXSENS/*RowiDV] . std(MAXSENS/^Row8DV 1 , tap; 
/•••• initialixing*^**''**''****/ 
£or(i=0;i<Ro%rtD;i++) ( 

mnCil-O.Oi 

8tdlil=0.0; 

} 



^•••••CAlculating the mean •••••*/ 

for »i=0 : i<ColsD; i**) 

f or ( j = 0 : j <RowbD ; j ) 
innlj]*-D(j] li) ; 
for ( j=0 ; j <RowsD; j**) 

ami j] /»ColsD; 

/•••• Calculating the standard deviation' / 
for (i"0; i<Col8D; i**) 

£or( j"0; j<RowaD; jf*) 

■ ,td(j]*-lD(jlti)-innlj])MD[j)[il-!nn(D]l; 

fort j=0; j<Row8D; j-^*) 

std(j}/><ColsD-l) ; 

/•••• Extending the range of the training matrix***/ 
for(i=0;i<ColsD;i+*) ( 

for( j»0: j<RowsD; j**) ( 

if (Dtjl [il<-fnn(j]) ( 

Dlj) [il— itdtjl^cx; 
Dttil [j3— 8td[jl*ex; 

) 

else( 

D(jl (i]*-8td[j3*cx: 
Dtlil [jl*-8td(jl*ex; 

) 

) 

) 

int lengthFile(FILE •£!)( 
char c; 
int i; 
i«0; 

while( (c«g«tc(fl) ) !«EOF) 

if(c»«'\n*) i**? 
£se«)c(£1.0L.O) ; 
return < i ) ; 

} 

/• 

int widthFile(FILE •£!)( 
char c; 

int i; 
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void Equaj.lySpaced( ) { 

int L,W,i. j ,k, count, nim. estop; 
float ex; 

InitialScreenl 1 ; 
t«xtLoo>t (YELLOW. BLUE) ; 
£aeek(Train.OL.O) ; 

90toxy (17,8) ; 

cprlntfC train using the full file? ly/n) : 
gotoxy(56, B) ; 
cvgetche ( ) ; 

i£(c»»'n' ) { 

gotoxy(17,8) ; 

cprintfC How many lines into the file: 

gotoxy(53. 8) ; 
cscanf <'%d-,tL) ; 

) 

elte( 

L»l©ngthrile (Train) ; 
f seek (Train. OL.O) ; 

) 

StOp«0 ; 

while (istopM 

gotoxy(14.8);/* */ 

cprintfC Eater nuacr .of training vectors (max 132): *); 
gotoxy ( 61 . 8) ; 
cicanf (•%d'.tnum) ; 
if(nun>132) 6top"l; 

} 

gotoxy (14, 8 ) ; 

cprintfC Getting vectors from training file.... "); 

count=floor( ( (float ) LI /( (float) nua) ) ; 

if(count<l)( count»l; L«n\2m; ) 

ColsDsnum; 

D-0; 

for ( iaO; i<L; i+*) { 

i£( (i%count)e«0) ( 

W>readinputSample (Train) ; 
RowaD«W; 

f or (k=0 ; )c<RowsD; k++) { 

D(kl [j]=Sanipie(k) ; 
Dt[ jl tk]=Sainple(k] ; 

) 

) 

) 

gotoxy(17.8) ;/• */ 

cprintf(' Extend training range (y/n) : •); 

gotoxy(56. 8) ; 

/•C8can£C%c*,fcc) ;*/ 

getche () ; 

cgetche ( ) ; 

if (c— *y ) ( 

gotoxy(17,14) ;/* */ 
cprintf ( • How many standard deviations from the mean: ■>; 
gotoxy(17. 15) ; 

cprintfC (typical value « 3) Enter value ■ •); 
gotoxy (5 3. 15 ) ; 
cscanf C%f ' .bex) ; 
gotoxy(17,17) ; 
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nuinPairs»2*W; 
didSwicch«l ; 
while (didSwitch) { 

numPaira-- ; 

didSwitch-0; 

forti=0:i<nuinPairs;i++) ( 

if (tiine_pcs(il>tiine_pcsli*ll ) ( 
tinp=tiinejpcs ( i] ; 
tiine_pta I i ) «tinie j)ts ( i+1 1 
tiine_pts ( i+1 ) =tinp; 
didSwiCchsl ; 



) 



) 



getting unicjrue time points ••••••/ 

gotoxy(17.8) ;/• . ^ 

cprintfC Removing redundant txme points. 

j»=0; 

forti«l:i<2*W;i*+) { 

if (time_pts(i-ll ! i=time jts [il ) { 

tiine_ptst j]-time_pts(il ; 

) 

} 

ColsD= j ; 
Row3D«W ; 
j«0; 

gctoxy(17,8) ;/* 

cprintf (• Creating training matrix 

for(i = 0; i<L;i+-«') ( 

w=readinputSainple (Train) ; 

if (i==time_pts[j] ) ( 

for (lc=0;k<RowsD;)c++) { 

Dtk] [j]»SainpleIk] ; 
Dt( j) ()c] -Sample [k] ; 



} 



) 



Extend training range (y/n) : 



) 

gotoxya7,8) ;/* 

cprintf ( • 

gotoxy (56 , 6 ) ; 

/•cscan£C%c'.fcc) ;*/ 

gecche ( ) ; 

c»getche ( ) ; 

if (c=='y) ( 
gotoxy (17 , 14) ; /• , 
cprintf (• How many standard deviations from the mean. 

gotoxy(17.15) ; , 
cprintf (■ (typical value = 3) Enter value • 
gotoxy (53 , 15) ; 
cscanf (•%£'. fc«x) ; 
gotoxy ( 17. 17 ) ; 

cprintf (■ Extending training matrix 

cxtendD(ex) ; 

} 



•); 



■1 
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/• KINIMAX TRAINING MODULE 



void MinMaxO { 



int L,W, j,i,)c,numPairs,didSwitch; 
int timej)ts(MAXSENS»2] ,c; 
/•int min_t:iine_pcs{MAXSE2JSl ;•/ 
floac minixiax(2) [MAXSENSl .trap, ex; 

InitialScreenO ; 
texcLook (YELLOW, BLUE) ; 
fa«ek (Train, 0L,0) ; 

gotoxy(17,e) ; 

cprintfC train using the full file? (y/n) : •); 
gotoxy (56 , 8) ; 
cs9ecche( ) ; 
/•cscanf ( ■%c" ,&c) ;*/ 

if(c««'n*)t 

gotoxy(17,8); 

CprintfC How many lines into the file: *) , 

gotoxy{53, 8) ; 
cscanf (•%(!', tL); 

) 

elseC 

L«lGngthFile (Train) ; 
£Boe)c (Train, 0L,0) ; 

} 

W=readinputSanTple (Train) ; 

gotoxy(17,8);/* •/ 
cprintfC Finding max and min for all sensors... '); 
for( j«0; j<W; { 

minmaxtO] [ j ) «Sainple [ j ] ; 

minroax(l] [ j ] ^Sample [ j 1 ; 

time_pts ( j]»0; 

tixne^pts [W+j] «0; 

} 

Finding time points for each sensor ••♦•••/ 
for (i=l; i<L; i++) { 

W«readinputSample (Train) ; 
for ( j=0; j<W; j**) { 

if (Sample(j)>«minroax[0] Ij] ) 

tifflej)t6[j)«i; 
if (Sanqjlel jl<«Tidnmax(ll (j] ) 
tinie_ptfl[W-l*j)«i; . 

} 

} 

f seek (Train, 0L,0) ; 

gotoxy(17,81 */ 
cprintfC Sorting time points "); 

/ / 

/••••••••sorting ••••♦♦•/ 
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Appendix B 

Computer software for producing an optimal uaining set derived by searching 
signal information during a training period to construct training vectors including 
highest and lowest points for signals under surveillance. 
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void lliplr (float a[), inc size) 

( 

int i , j ; 
float rtmp; 

3-floor isiie/2) ; 
forU-0;i<j;i++) 
{ temp«a [i] ; 

a[i]-a[6ize-i-ll ; 

a laize-i-l) -temp; 

) 

) 

void dertfloat y[]. float dl£fy(] , int size) 

{ 

int i; 
float diff; 

for(i-0;i<flize-l;i**) 

diffytx] - y(i*l]-y[i] ; 

) 

void rootstfloat dl£fxc(], iat size, float Toot) 

( 

int 

float term; 

if (diffxc(03*diffxc(size-l] > 0) 

{ printf("The cross^correlation is not unimodal . \n") ; 
exit(0}; 

) 

else 

{ (•root)«0.0; 
for (i«0; i<si2e-l ; i*4-) 
{ termsi; 

f or ( j«0; j<8izc-l; j^-f) 
if ((i-j)>0) 

term-term* ( (-dif fxc ( j] I / (dif fxc [i] -dif fxc t j] 1 ) ; 
(*root) - ('root) term; 

1 

) 

} 



void pm info (char •pgin_name) 

( 

printf ("\n%8%B%s\n", 

"Usage: pgm^name. " inlllel infile2 in£ile3 outfileM ; 
printf("\n-)? 

printf("Note that infilel & infile2 are reference and input signal"); 
printfC* respectively. in£ilc3 should contain the coefficients of); 
printfC the LPF, and outfile receives the phase •equalized input signal. Vn") ; 

) 
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ph«ti.l_xc[0) ; 

j-0; " . . 

for(i.O;i<fil.xcii«:i**' 

if (fil.xcCi)>ph) 

{ ph«fil_xc(il ; 

j-i; 



) 



priiit£{-The Bhift is %d\n-.ph): 

/• Find the shift, U5xng derivative technique . 

if (method— 1) 
( der(£il xcdifffil.xc.fil^xcsixe). 

rootB (dif f f il^xc. £ll.xc8ixe-i . troot) ; 

ph-root- (fil^xcBitetl) /2; 

printffThe shift is %f\n-.ph); 



) 

/• phB.e-equaIixe the input (f2) relative to reference (fl) 
If ph>0.^ leads fl; if ph<0, t2 lags fl. */ 

if lph>0) 
{ forU-ph;i<N;i+*» 
f2(i-phl - f2li]: 
forti-N-ph;i<N;i++) 
f2lil-0; 

) 



if (ph<0) 
{ ph • -ph; 

for(i»N-l;i>-ph;i--) 

f2(i} - f2(i-pb] ; 
for (i-0;i<ph;i*+) 
f2li]-0; 

) 

for(i-0;icK;i+*) 

fprintf (outf ile . "%£\n" . f 2 til ) ; 



f close (infilcl) ; 
f close (infile2) ; 
fclose(in£ile3) ; 
f close (cutfile) ; 

) 



void conv( float fll. float gtl". float ell. 
{ 

int m.)c; 

for ( k-0 ; k<8irc ; k** ) 
{ cW-0; 

f or i tn« 0 ; m< -Jc J tn+ ♦ ) 
cM « c(kl+f ltnl*g(k-m) ; 

) 

) 
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/• open the input and output data files, •/ 

outf ile-£openiargv{4l , "WM ; 

f or li-0;i<xcsi2e; 
( fllil-O; 
f2[i].0; 

) 

for (i«0;i<f il_xcsi2e; i**) 
{ xcUl-0; 
b(il-0; 

fil_XC(x]-0; 

} 



if ( (infilel«fopen{argv(l) , -r") ) KULL) 
{ printf ("There is no data file %8t\n", argv(l]); 
exittO) ; 

) 

else 

{ i»0; 

while { (f scanf (inf ilel , . fi.fl [i]) ) ! =.EOr) 

) 



if ( {infile2«fopentargv(21 , -r-) ) NULL) 
{ printf (**There is no data file %9!\n", argv[2]); 
exitlO) ; 

) 

el»e 

{ i-0; 

while ( (f scanf (inf ile2 , nf - , tf 2 [i] ) ) ! -EOF) 

} 

if ( (infile3-fopentargvt3l , "r") ) «- NULL) 
{ printf ("There is no data file %Bi\nV argv{3]); 
cxittOJ ; 

} 

else 

{ i-0; 

while { (f scanf (inf ile3 , " Vf , tb [il ) 1 i -EOF) 
i+4-; 

} 

/• Calculate the cross -correlation of the input and the */ 
/♦ reference signals. ♦/ 

flipir(fl,N) ; 

convtf 1. f2.xc,xc8ize) ; 

/•To pass xc through a LPF */ 
conv(b.xc,fil_xc. fil^xcsixe) ; 

Find the shift, using direct *max. correlation tecnnique. 
If (metnod«-2 ) 
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